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ABSTRACT 

We show that the differential rotation profile of the solar convection zone, 
apart from inner and outer boundary layers, can be reproduced with great accu- 
racy if the isorotation contours correspond to characteristics of the thermal wind 
equation. This requires that there be a formal quantitative relationship involving 
the entropy and the angular velocity. Earlier work has suggested that this could 
arise from magnetohydrodynamic stability constraints; here we argue that purely 
hydrodynamical processes could also lead to such a result. Of special importance 
to the hydrodynamical solution is the fact that the thermal wind equation is 
insensitive to radial entropy gradients. This allows a much more general class 
of solutions to fit the solar isorotation contours, beyond just those in which the 
entropy itself must be a function of the angular velocity. In particular, for this 
expanded class, the thermal wind solution of the solar rotation profile remains 
valid even when large radial entropy gradients are present. A clear and explicit 
example of this class of solution appears to be present in published numerical 
simulations of the solar convective zone. Though hydrodynamical in character, 
the theory is not sensitive to the presence of weak magnetic fields. Thus, the 
identification of solar isorotation contours with the characteristics of the thermal 
wind equation appears to be robust, accommodating, but by no means requiring, 
magnetic field dynamics. 

Subject headings: convection — instabilities — Sun: rotation — Sun: helioseis- 
mology 
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1. Introduction 

In a recent paper, Balbus (2009, hereafter B09) argued that the shape of the isorotation 
contours in the solar convection zone (SCZ) may be understood as a consequence of a dom- 
inant thermal wind balance in the vorticity equation, together with the near coincidence of 
surfaces of constant specific entropy S (isentropes) and surfaces of constant angular velocity 
f2 (isotachs). With S a function of f2, the thermal wind equation (TWE) may be solved 
analytically and exactly; the isorotation contours then explicitly correspond to the charac- 
teristics of the TWE. B09 found that the agreement between the characteristics and the 
observed Q profile is very good, even under the simplest of assumptions. In this paper, we 
take a more systematic approach, and find that the agreement is truly remarkable (cf. figure 
[1]). The quality of the fit is strong evidence not only that thermal wind balance holds in the 
bulk of the solar convection zone (as noted previously in many numerical simulations), but 
that there is also some functional relationship between entropy and angular rotation that 
needs to be understood. 

The hope, of course, is that this is telling us something important about the physics of 
the solar convection zone. B09 puts forth the case that the SCZ is in a state of marginal 
dynamic instability not only against convective instabilities, but against magneto-baroclinic 
instabilities. To understand why the distinction is important, recall that a weak magnetic 
field in a differentially rotating gas is a catalyst for destabilization. This effect is well 
known to the accretion disk community in the form of the magnet orot at ional instability 
(MRI), which renders fully turbulent what would otherwise be hydrodynamically stable 
Keplerian disks. B09 argues that in its magnetobaroclinic guise, this weak field instability 
can destabilize what would be hydrodynamically stable stratified configurations of the SCZ, 
driving the system to a condition of marginal instability. Such a state is reached when 
constant entropy surfaces and constant angular velocity surfaces nearly coincide, leading to 
an5 = f(Q 2 ) relationship. 

In this paper, we examine another possible explanation for the striking fit evidenced 
by figure [1] in §2. We are motivated, in part, by the fact that some purely hydrodynamic 
calculations are able to reproduce rather well certain noncylindrical features of the solar 
rotation profile (Miesch, Brun, & Toomre 2006, hereafter MBT06; Miesch et al. 2008; Miesch 
& Toomre 2009). Those runs show that, provided there are sufficiently large latitudinal 
entropy gradients at the base of the SCZ, a very plausible fit to the solar isorotational 
contours can be found throughout the bulk of this region. In this current paper, we will 
argue that there may well be a generic hydrodynamical mechanism that explains both the 
observed helioseismology results and the numerical simulations. 

An outline of our paper is as follows. In §2, we present an improved fit between an 
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explicit solution of the thermal wind equation and the helioseismology data. We also in- 
troduce, following MBT06, the concept of "residual entropy": the average entropy profile 
remaining after an underlying radial profile has been removed. Arguments are presented 
which suggest that convecting fluid elements will tend to move in surfaces of constant resid- 
ual entropy. In §3, we further argue that the same fluid elements will also tend to move in 
surfaces of constant angular velocity, and that these surfaces must therefore coincide with 
those of residual entropy. The coincidence of surfaces of residual entropy and angular veloc- 
ity is necessary in our approach to obtain a solution of the thermal wind equation. In §4 we 
conclude with an evaluation and discussion of the relative merits of the hydrodynamical and 
magnetohydrodynamical approaches. 

2. A simple model for the solar isorotation contours 
2.1. Preliminaries 

2.1.1. Coordinates and notation 

Let (R, <p, z) be a standard cylindrical coordinate system, and (r, 6, 0) a standard spher- 
ical coordinate system. We consider an equilibrium flow in a state of azimuthal rotation in 
which the angular velocity Q is assumed to be independent of 0, but may depend upon R 
and z. The background entropy profile S is also a function of R and z. Our notation for 
the other fluid variables is standard: v is the velocity, P is the gas pressure, and p is the 
mass density. Unless otherwise stated, all thermodynamic variables are understood to be 4> 
independent azimuthal averages. The velocity v will in general contain (convective) fluctu- 
ating components; any azimuthal averaging will be treated explicitly. The SCZ gravitational 
field g is to a good approximation GM Q /r 2 , where GM & is the product of the Newtonian 
gravitational constant and a solar mass. Finally, we define a dimensionless entropy function 
a: 

a = \nPp-\ (1) 

where 7 is the adiabatic index. The thermodynamic entropy is then given by S = Cpa, 
where Cp is the specific heat at constant pressure. 

2.1.2. Review of the TWE solution 

Our starting point is the partial differential equation arising from the assumption of 
a dominant thermal wind balance in the time averaged (and thus azimuthally averaged) 
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vorticity equation (e.g. Kitchatinov & Riidiger 1995, Thompson et al. 2003, Miesch 2005, 
B09): 

an 2 g da 

oz jr o9 

If, for some reason, there is a functional relationship of the form a = f{Q 2 ) (the sign of Q 
presumably does not matter), then the TWE may be written in spherical coordinates as 

dtt 2 ( gf tan0\ dQ 2 

1 ] 0, (3) 
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where /' = df/dfl 2 . Notice, however, that exactly the same equation obtains if the functional 
relationship between o and Q takes the more general form 

a' = a-a r = f(Q 2 ), (4) 

where a' will be referred to as residual entropy, and a T is any function of r alone. We will 
make important use or this "gauge freedom" later in this paper. 

The solution of equation ([3]) is that Q 2 is constant along the characteristic contours 
(B09): 

R 2 = r 2 sin 2 9 — A , (5) 

r 

where A is a constant of integration, and 

7 

If we denote by the subscript '0' the fiducial starting point of the characteristic (at which Q 
is specified), then the contour takes the form 

R 2 = r 2 sin 2 9 = rl sin 2 9 - B ( - - -) (7) 

While B is a constant along a given characteristic, its value can change from one char- 
acteristic to another depending upon the form of the function /'. B09 showed, however, that 
the simplest parameterizations of /' {e.g., f — a global constant or /' = a linear function of 
sin 2 6'o) already give excellent qualitative fits to the helioseismology data. In fact, it is not 
difficult to do even better. 



Figure 



1] shows a comparison between isorotational contours (in black) from recent 



GONG datcQ and the characteristics of the thermal wind equation (overlayed in white). 



1 We thank R. Howe for providing us with these beautiful results. 
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Here, we have carried out a more systematic approach than that of B09: we have integrated 
each characteristic from an interior starting point of 0.9 solar radii (as opposed to the less 
accurate surface fitting in B09), and chosen the constant B to match the initial slope of 
the data at the fitting radius. The new result is remarkably accurate, not just in qualitative 
form, but in quantitative detail as well. Expected departures of the black and white contours 
in the surface layers and tachocline are plainly visible, yet in the bulk of the SCZ the match 
is excellent. 

Note the equator and poles, where the curvature of the contours is, perhaps surprisingly, 
reproduced with great precision]^ While this is evidence of the validity of thermal wind 
balance, the implied functional relationship between Q and a (or a') is also at least partly a 
consequence of mathematical symmetry. In the polar neighborhoods R gradients are small, 
and near the equator z gradients are small, so that in both regions azimuthally averaged flow 
quantities will be a function predominantly of spherical r. These quantities must therefore 
be functionally related. This functional relationship evidently extends seamlessly into the 
regions immediately adjacent to the symmetry points, in which the isocontours are not 
predominantly radial. 

2.2. A functional relation between residual entropy and angular velocity 

2.2.1. Theoretical considerations 

One possible explanation for the striking resemblance of the analytic isorotation contours 
to the solar data was advanced in B09 by invoking a weak magnetic field. In the presence 
of such a field, disturbances associated with magnetobaroclinic modes lead to a condition of 
marginal instability if o = f(Q 2 ). This is a dynamical connection between a and Q, indeed 
an MHD connection. 

We will review the merits and shortcomings of this scenario in more detail in §4. Here, 
we explore a more general connection between a and Q relying neither upon magnetic fields, 
nor even upon dynamical constraints. The idea relies instead on the gauge freedom expressed 
in equation (TJJ. 

Consider a fiducial, nonrotating, spherical, convecting star that maintains a well-defined, 
long-term average radial entropy profile a r . Convection continously mixes high and low en- 



2 The inversion process by which the rotation contours are determined from helioseismology data is, it 
should be noted, less reliable near the poles. 
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Fig. 1. — Constant Q contours ([7]) of the thermal wind equation ([2]) (white curves) plotted on 
top of (black) isorotation contours from helioseismology data (GONG results courtesy of R. 
Howe). Blue contour is the bottom edge of the convective zone. Scale is in solar radii. Away 
from the tachocline and outer surface layers, the match is excellent. See text for further 
details. 
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tropy fluid elements, yet a r remains fixed. The reason for this is not mysterious. The system 
is heated from below, and this input tries to drive the gas into a state with a considerably 
steeper entropy profile. Convection counters this, allowing a balance to be struck between 
the external thermal driving and the heat transport by turbulent fluid elements. In particu- 
lar, convective displacements do not need to be within constant entropy surfaces to preserve 
the long-term a r profile. This profile is already being maintained by the external heating. 

If we now introduce a small amount of rotation into the problem, we expect Coriolis 
forces to produce a more complex, but still well-defined, long-term entropy profile, a(r,9). 
As before, convection will try its very best to change this profile. However, in contrast to 
the spherical case of the previous paragraph, it is very far from clear that radial heating 
from below can maintain a steady, nonspherical entropy profile. On the other hand, if the 
convective velocities v on average satisfy the constraint 



then the entropy profile certainly can be maintained. In this case, the long term maintenance 
of the spherically averaged component of the entropy o r is assured by heating from below, 
while the residual entropy a' is preserved because convective displacements will tend on 
average to occur in surfaces in which this quantity is constant. Thus, the long-term average 
entropy profile a(r, 9) remains intact. Convecting elements move not in constant entropy 
surfaces (which would preclude heat transport!), but in surfaces of constant residual entropy 



Some care is needed for the proper interpretation of equation (JSJ) . The poloidal compo- 
nents of v (hereafter v p ) have a long term average value much smaller than a typical fluctu- 
ation amplitude \v p \. (This systematic velocity manifests itself principally as the meridional 
flow that has been measured at, and near, the surface.) A time average of equation (jSJ) 
would produce contributions of second order in the fluctuations, including heretofore ne- 
glected a' fluctuations. On the other hand, the fact that the "up-down" convective velocity 
fluctuations nearly cancel means that for a given convective cell, there is very little to dis- 
tinguish these two directions. (Because of rotation, "upward" and "downward" need not 
be precisely radial, of course.) These first order velocity fluctuations in essence establish a 
well-defined axis relative to which Vcr' is orthogonal; the positive and negative sense of the 
axis is unimportant. A more precise formulation of equation (jSJ) is therefore 



where the angle brackets ( ) denote time averaging. This is also a more robust formulation, 
since the quantity on the left is now dominated by the first order velocity contribution v p , 
and fluctuations in a' do not enter. 



v . Vff' = 





(9) 
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More physically, simulations and laboratory experiments suggest that convective trans- 
port is dominated by long-lived coherent structures. The updrafts and downward plumes 
characteristic of these structures are clearly associated with dominant "first order" velocity 
terms, and it is these structures that would exist within and maintain constant a' surfaces. 
If it could now be further argued that convecting fluid elements (or coherent structures) 
also lie within constant Q surfaces, we would arrive at an understanding of why there is a 
functional relation between Q and a': if elements are simultaneously moving in r6 surfaces of 
constant Q and a', these must be the same surfaces. In other words, the two quantities are 
functionally related. In the next section, we will see that there is indeed a simple argument 
leading to the conclusion that displaced fluid elements will inevitably move in constant Q 
surfaces. 

Finally, we reiterate that although the original B09 derivation of the isorotation contours 
was based on the assumption that a = /(fi 2 ), exactly the same contours emerge if the func- 
tional relation is a' = /(fi 2 ). This very simple gauge invariance of the thermal wind equation 
allows for a considerably wider, more general class of thermal wind solutions, eliminating 
the constraint that constant (total) entropy and angular velocity surfaces coincide. 

2.2.2. Comparison with simulations 

Support for the line of argument advanced in the previous section can be found in 
the (purely hydrodynamical) MBT06 simulations. In these runs, the residual entropy a' is 
calculated relative to a background radial profile, corresponding to our a r . In figure [2], 
we show time averaged contours of the angular velocity profile Q from the MBT06 run 
designated AB3, of the residual entropy, and of the total entropy a. The agreement between 
the Q and a' contours is apparent, as is the gross disagreement between fl and a. Evidently, 
there is a functional relation between a' and Q, whereas there is no such relation between 
f2 and a. This is in good agreement with the above discussion. Notice as well the small 
gradient of a' near the equator. This is as expected if a r represents an angular average of a, 
which would be dominated by 6 values near ir/2. 

3. Nonaxisymmetric disturbances in a differentially rotating medium 

3.1. Comoving coordinates and wavenumbers 

Consider the behavior of convective disturbances in a differentially rotating gas Q(R, z). 
The disturbances are embedded in this shearing medium, and since they are not propagated 
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Fig. 2. — Isorotation contours (a), contours of constant residual entropy (b), and contours 
of total entropy (c), from the run AB3 of MBT06. (In [c], a constant background entropy 
has been subtracted out.) There is generally very good agreement between the angular ve- 
locity and residual entropy contours, but not between the angular velocity and total entropy 
contours. (Figures courtesy of M. Miesch.) 
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as waves, we may think of them as local in character. But "local" refers to the shearing 
gas, not to an absolute space-time grid, so it is particularly useful to express flow quantities 
in Lagrangian coordinates tied to the shear flow. We designate these coordinates as the 
"primed" system. This transformation is given by 

# = R, ft = <j> - n(R, z)t, z' = z, t' = t. (10) 

For our present purposes, we need consider only the equation of mass conservation: 

V-(pv) = 0, (11) 

where v is understood to be the velocity relative to the differential rotation RVl. The 
gradients are of course the standard Eulerian derivatives taken at fixed time. We have not 
written the explicit p time derivative, since we are ultimately interested in long term time 
averaging and the term will disappear. Apart from this, the equation is exact; the mass flux 
is taken to have a full space time dependence. 

Transforming to our comoving system, 

ao r) 

(12) 
(13) 
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dQ d 
t dR dft' 
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dtt d 
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dz 





and the derivative is unchanged. More compactly, 

v = V-(tvn)^, (14) 

a relation that holds for all three components. 

The embedded mass flux field pv may be expressed as a superposition of Fourier com- 
ponents of the form 

exp [i (k' R R' + m(f)' + k' z z' - cut')] (15) 

in which all components of the k' wave vector are constants. The above coordinate trans- 
formation implies that in the mass conservation equation, the standard Eulerian spatial 
derivatives with respect to R and z are replaced respectively by ik R (t) and ik z (t), where 

k R (t) = k' R -mt—, (16) 
dQ 

k z (t) = k' z -mt— (17) 

Henceforth, the time dependence of k R and k z will be understood. A somewhat similar 
formalism has been developed for embedded disturbances in shear turbulence, where it is 
known as rapid distortion theory (Townsend 1980). 
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3.2. Mass flux morphology 

The Fourier series for the mass flux is: 

pv = ^ v( k ', u ) ex P i k> R r ' + m $ + K z> ~ , ( 18 ) 

where fi(k', uj) is the amplitude of the indicated wavenumber and frequency. Mass conser- 
vation leads immediately to 

= V-(pv) = k-n{k', u) exp [i (k' R r + m(j)' + k' z z - ut')] . (19) 

Each Fourier term is linearly independent of the others, hence for all k', u, 

(k' -mtVtt)-n(k',uj) = (20) 

For nonvanishing m, at sufficiently large t, 

/i(fe»-vn = o. (21) 

The mass flux is a superposition of the /j, vectors. Therefore, if the flux is not dominated by 
an m = component, it will tend toward surfaces of constant Q. Note that this statement 
is restricted neither to linearized displacements, nor to the WKB limit. Finally, it has been 
noted that this tendency for turbulent convective structures to align with constant Q surfaces 
may also be supported on a dynamical basis via the effects of the Coriolis force on plumes 
(Brummell, Hurlburt, & Toomre 1996). 

Alignment of embedded structures with constant Q surfaces is of course seen regularly 
in accretion disk simulations, where the shear is large and depends only upon R. Disks 
features often appear nearly axisymmetric because of the effect of the shear. In the Sun, the 
R and z dependence of Q leads to a more structured response. This assumes, of course, that 
the rotational shear is strong enough to interact with the convective cells over the course of 
their lifetime. This is a reasonable assumption for the Sun, especially if the convection is 
dominated by long-lived coherent stuctures. 

This, in principle, supplies the missing link of the argument of §2.2.1: the course of 
a displaced mass element will tend on average to follow a constant Q surface, essentially 
because differential rotation wraps the flow into sheets of constant Q. Since fluid elements 
move at once in surfaces of constant fl and constant a', these surfaces must in general 
coincid^E 



3 The degenerate case in which the surfaces do not coincide would correspond to trivial circular azimuthal 
orbits of the fluid element. 



- 12 - 



Notice how the complete set of three hydrodynamical equations — mass conservation, 
entropy conservation, and the equation of motion — has been incorporated in our solution. 
From mass conservation, we infer that fluid elements move in constant Q surfaces. From 
entropy conservation, we infer that fluid elements move in surfaces of constant a'. From 
these two inferences, it follows that constant Q and constant a' surfaces coincide, and this, 
expressed as a functional dependence, is used in the equation of motion (thermal wind 
balance) to construct the isorotational contours of figure [1]. 



3.3. Heat transport 

The convective heat transport is 

Q = pv5w, (22) 

where Sw represents the difference between the specific enthalpy of a convecting element and 
its surroundings. The reasoning we used in the last section shows that 

Q ■ Vn -> 0, (23) 

as time progresses. In other words, the heat flux, like the mass flux, is predominantly in 
surfaces of constant f2. 

This offers a very plausible physical basis for what has long been regarded as a puzzling 
issue: why are the solar rotation contours quasi-radial? The answer we suggest is that the 
rotation contours are also the natural conduits for convective heat transport, and their near 
radial character simply reflects this property. The contours represent a compromise between 
strong radial thermal driving and the dynamical exigencies embodied in the TWE. Thus, 
they are nearly, but not perfectly, radial throughout the bulk of the SCZ. Near the equator 
and poles, the convective heat flow certainly does not follow isorotation contours, but as 
has already been noted in §2.1.1, these are precisely the regions in which the functional 
relationship between a' and Q is maintained by mathematical symmetry. 

It is interesting to speculate on how the isorotation contours would change in rapidly 
rotating stars. For example, when a solar-type star first arrives on the main sequence it spins 
much more rapidly than the Sun does now. The convection pattern is likely to be dominated 
by Coriolis effects, expressed through the Taylor-Proudman constraint (Busse & Simitev 
2007). Numerical models confirm that convection takes the form of elongated "banana cells" 
at low latitudes, outside the tangent cylinder that encloses the radiative zone (Ballot, Brun 
& Turck-Chieze 2007; Brown et al. 2008). In these rapidly rotating stars, the isotachs tend 
to be roughly cylindrical and there are pole-equator differences in entropy and temperature 
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that are significantly larger than those inferred for the Sun. In our model, the 1/Q 2 scaling 
of f'(Q 2 ) suggests a much smaller value for the B parameter in rapid rotators, which does 
indeed lead to more cylindrical isotachs. 

4. Concluding discussion 

What is novel about our approach is not the use of the thermal wind equation per se, 
which has provided important guidance to workers in this field for many years now. It is the 
suggestion that there is a functional relationship between the entropy and angular velocity, 
for it is this that allows an analytic deduction of the shape of the Sun's isorotation contours 
to be made directly from the TWE. 

As a formal matter, there is some freedom in choosing the precise form of the entropy- 
angular velocity relation, a type of gauge invariance. This is because the TWE involves only 
the 9 derivative of the entropy and the z derivative of the angular velocity. Different physical 
mechanisms for coupling a and Q in principle entail a different choice of gauge. 

In B09, an MHD process was invoked to couple a and Q. That paper noted that 
baroclinic modes in a weakly magnetized gas become marginally unstable if constant o and 
constant Q surfaces coincide, i.e. a = /(fi 2 ). In that approach, rotation and stratification 
are set on an equal footing in determining the dynamical stability of the Sun's outer layers. 
Moreover, in the MHD model, the radial entropy gradient is restricted to be significantly 
less than the 6 gradient. This can be readily deduced from figure [1], since in the a = /(fi 2 ) 
interpretation, the GONG data curves are isentropes as well as isorotation contours. Whether 
the actual radial entropy gradient is constrained in this way depends upon the scale of the 
effective mixing length. If this length is only a small fraction of the radial extent of the 
convection zone, the radial entropy gradient could significantly exceed the 9 gradient. 

Finally, in the MHD model, the goodness of the fit between the theoretical and observed 
isorotation contours depends upon a very tight functional relationship between o and Q. 
But were the a and Q surfaces to coincide exactly, convection would be completely cut- 
off. Conversely, as we have noted above, if very vigorous convection is required, it becomes 
untenable to maintain a a = f{Q 2 ) relation. Is the departure between Vcr and Vfi large 
enough to allow for vigorous convection, yet small enough to ensure that o = /(fi 2 ) remains 
an excellent approximation in the TWE? 

The purely hydrodynamical explanation for the solar isorotation contours avoids many 
of these potential difficulties. The restriction on the size of the radial entropy gradient is 
lifted by seeking a functional relationship of the form of equation As the residual entropy 
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o' can be significantly smaller than a, the effects of a potentially large radial entropy gradient 
can be eliminated from the rotational dynamics of the problem. The physical basis for this 
functional form is related to the fact that convection should occur on average in surfaces of 
constant a', so as to maintain a well-defined time averaged a (total entropy) profile. The 
kinematics of differential rotation and simple mass conservation ultimately confine mass 
motions to surfaces of constant Q. Together with equation ([8]), this ensures a relationship 
of the form a' = /(fi 2 ), and even suggests why the contours are quasi-radial (cf. eq. [23]). 
Without the need to restrict da/dr, convection can be as vigorous as needed (provided that 
it does not dominate the differential rotation), yet not directly affect the TWE solution for 
the isorotation contours. 

Finally, the "hydrodynamical" argument can easily accommodate MHD processes. Nei- 
ther mass conservation nor entropy conservation, the two physical processes central to our 
argument, change in the presence of a magnetic field. The dynamical details, in particular 
the rotational couplings, may well be altered, but the functional relation between a' and Q 
should remain robust. We therefore suggest that the simple physical arguments leading to 
equation (j4j) may prove to be a suitable basis for understanding much of the dynamics that 
gives rise to the solar isorotation contours. 

The theory outlined in this paper is meant to be a viable alternative to an MHD model, 
and it is as yet far from complete. Our presentation falls well short of the standards of 
mathematical proof. We have relied on "tendencies" for the fluid to move in surfaces of 
constant Q and o' to be strong enough to firmly establish the relationship expressed com- 
pactly by equation (J4]). But figure [1] speaks for itself. To understand at a deeper level why 
this approach seems to work so well, a more detailed elucidation of the dynamics is clearly 
desirable. 

To conclude, in this paper we have advanced theoretical arguments to explain the strik- 
ing correspondence between isentropes and isotachs in azimuthally averaged models of the 
solar convection zone. These arguments need to be backed up by directed and detailed com- 
putations. To this end, we are currently developing an anelastic code that can be used to 
demonstrate explicitly that the relationship of equation (J3J) is valid for a variety of relevant 
configurations. Our principal aim will be to explore the connection between Q and a (or a') 
as the overall rotation rate is varied. The effects on this relationship of varying the thermal 
boundary conditions, in order to represent the influence of the solar tachocline (Miesch & 
Toomre 2009), are also of great interest. 
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